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ABSTRACT 

The collapse of an isolated, uniform and spherical cloud of self-gravitating particles 
represents a paradigmatic example of a relaxation process leading to the formation of 
a quasi-stationary state in virial equilibrium. We consider several N-body simulations 
of such a system, with the initial velocity dispersion as a free parameter. We show 
that there is a clear difference between structures formed when the initial virial ratio 
is &o = 2Kq/Wo < b^ w —1/2 and bo > fog. These two sets of initial conditions 
give rise respectively to a mild and violent relaxation occurring in about the same 
time scale: however in the latter case the system contracts by a large factor, while in 
the former it approximately maintains its original size. Correspondingly the resulting 
quasi equilibrium state is characterized by a density profile decaying at large enough 
distances as r~ 4 or with a sharp cut-off. The case bo < fog can be well described by 
the Lynden-Bell theory of collisionless relaxation considering the system confined in 
a box. On the other hand the relevant feature for bo > fog is the ejection of particles 
and energy, which is not captured by such a theoretical approach: for this case we 
introduce a simple physical model to explain the formation of the power-law density 
profile. This model shows that the behavior n(r) ~ r~ 4 is the typical density profile 
that is obtained when the initial conditions are cold enough that mass and energy 
ejection occurs. In addition, we clarify the origin of the critical value of the initial 
virial ratio fog. 
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1 INTRODUCTION 

The evolution of a system of massive particles interact- 
ing solely by Newtonian gravity is a paradigmatic problem 
for astrophysics, cosmology and statistical physics. The un- 
derlying open question concerns the relaxation mechanism 
that drives the system to form structures which seem to 
be in a sort of equilibrium, as for instance different kind 
of astrophysical objects such as globular clusters, galax- 
i es, and galaxy clusters (JLvnden-B cll 1967; Pad manabhanl 



1990; Binncv & Trcmainc 1994; Saslaw 2000: Hcggies 2002 



Aarsethl 120031 ). In a galaxy th e two body relaxation tim e 



is of order r2 ~ 10 17 years (jBinnev fc Tremaind 1 19941 ). 
and is much longer than the age of the universe (i.e., 
~ 10 10 years): for this reason these objects are not in 
thermal equilibrium. However, they prese nt common fea- 
tures as the luminosity profiles (see e.g., Ide Vaucouleursl 
{ 19481 ) : lBinnev fc Merrifieldl (|l998l )). Much theoretical work 
has been devoted to study the dynamical model to character- 
ize such profiles and despite the numerical simulations have 



shown that structures formed in some cases are compati- 
ble with observations, the physical origin of these profiles 
has not been yet clarified from a theoretical point of view. 
Namely, the problem still remains to explain how to form 
the shape of density profiles and of velocity distributions of 
stellar structures like elliptical galaxies and globular clus- 
ters that are generally characterized by a dense central core 
and a dilute halo — where the halo is often featured by a 
power-law decay of the radial d ensity (jBinnev fc Tremaind 
ll994l : lBinnev fe Merrifi"eTdlll998l ). 

In cosmology one faces a different but somewhat related 
problem. Since more than a decade it has been realized that 
a major issue about gravitational clustering dynamics con- 
cerns the formation of the so-called halo-structures, which 
are considered the primary building blocks in terms of which 
the non-linear struc tures observed in cosm ological simula- 
tions are described (jCoorav fc Shethl 120021 ). These are ap- 
proximately spherical symmetric structures, but sometimes 
with complex substructures, and with a density profile that 
that has almost universal statistical features and unknown 
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dynamical origin. Density profiles of dark matter halos have 
become one of the most challenging issues for our under- 
standing of cold dark matter structure formation. Numer- 
ical simulations provide evidences of steep central density 
cusps with power law slopes p ~ r~ , with j3 ~ 1 at small 
scales and ~ 3 a t larg e ones l|Navarro et alJll99a. 11997c 



Moore et al. 
20051: 



KK- 



Navarro et al 



iDiemand etal.| J2004 iReed etal.1 
20041 ; iMerritt et alj|200c^ . Recently 



200 



Graham et al] ([20061 ) showed that in simulated dark matter 
models, at large enough scales, slopes of j3 « —4 might be 
permitted. Several attempts have been made f or an analyt- 
ical d eri vation of the d e nsity p rofile (see, e.g.jBertschingerl 
dl985h; ISver fc White! Jl99ST>; I Subramanian et al.l tool) ; 
Iffiotelisl | |2002J ); iManrique etal.1 (|2002l ); iDekel etal.1 Ift&M ) 
and references therein), and none seem to present a clear 
and simple explanation for the findings of N-body codes. 

The question of the nature of the equilibrium properties 
of these core-halo structures is thus relevant both in astro- 
physics and cosmology and thus one would like to develop 
a statistical mechanics approach to describe these systems. 
However, one must consider that, from the point of view of 
statistical physics, self-gravitating systems present funda- 
mental problems, that are also common to other long-range 
interacting systems. Indeed, it is well known since the pio- 
neering works of Boltzmann and Gibbs, that systems with 
a pair potential decaying with an exponent smaller than 
that of the embedding space, present several fundamen- 
tal problems that prevent the use of equilibrium statisti- 
cal mechanics: thermodynamic equilibrium is never reached 
and the laws of equilibrium thermodynamics do not ap- 
ply ( Padmanabhanll 19901 ; iDauxois et al.ll2003l ; ICampa et al.l 
12008 ). Rather these systems reach, driven by a mean-field 
collisionless relaxation dynamics, quasi-equilibrium configu- 
rations, or quasi-stationary state (QSS ), whose lifetime di 



erge with the number of particles N (IDauxois et al] 



Campa et all |2008t lYamaeuchil 120081; iGabrielli et al] 



2003 



2010 



Joyce fc Worrakitpoonponll2012l ; IWorrakitpoonpon fc Joyed 



20121 ). The formation of QSS is at present one of the most 
living subjects in non-equilibrium statistical physics and a 
general theoretical framework is still lacking: it is thus nec- 
essary to consider toy models and/or relatively simple sys- 
tems that can be studied through numerical well-controlled 
experiments. 

In order to understand the formation of a core-halo 
structure, a paradigmatic example is represented by the col- 
lapse of a spherical, isolated and uniform cloud of N ran- 
domly placed particles with mass density po interacting only 
by Newtonian gravity. This s ystem has been considered since 
the early numerical studies (|Henonll 19641 ; Ivan Albadalll982l ) 
when it was realized that it relaxes violently, in a typical 
time scale td = \/37r/(32Gpo), to produce a virialized state. 
Such a time scale is much shorter tha n the two-body col- 
lisional time scale T2 « N/log(N)rD (JBinnev fc Tremainj 
Il994l ; ISaslawl 120001 ) and for this reason in the time range 
td < t < T2 the system relaxes into a QSS in virial 
equilibrium. Then, because of two-body collisions parti- 
cles can gain some kinetic energy and evaporate from 
the system: on a time scale of the order of ti the system 
changes shape because of particles evaporation. Simple con- 
siderations bas e d on the microcanonical entropy (see e.g. 
IPadmanabhanl |l99y)) imply that at asymptotically long 
times, and for a purely Newtonian potential, the particles 



will tend to a configuration in which there is a single pair 
of particles with arbitrarily small separation, and the rest of 
the mass is in an ever hotter g as of free particles so to con- 



serve the total energy (see e.g. lAarsethl (|l974T ); Ijovce et al] 
J2BQ3)). 

The underlying physical process in the formation of 
core-halo structures in the cosmological context is thought 
to be similar to the collective relaxation of such a finite 
and i solated self-gravitating particle system. iLvnden-Belll 
l| 19671 ). who named the collective relaxation process as "vi- 
olent relaxation" made a theoretical attempt to explain the 
gravitational collapse by approximating the temporal evo- 
lution as governed by the collisionless Vlasov equation and 
thus neglecting binary collisions. By introducing a coarse- 
graining in phase space the equilibrium state is postulated 
to be the one that maximizes the entropy computed by 
counting all the possible micro-states compatible with the 
Vlasov-Poisson conservations laws. In this context, differ- 
ently to ordinary thermodynamic equilibrium states, the sta- 
tistical properties of the QSS depend on initial conditions. 
The predictions of the Lynden-Bell approach were however 
shown to be at odds with the r esults of numerical experi- 
ments |Arad fc Johansson! I2005T ). The failure of the theory 
was attributed to the fact that the violent relaxation occurs 
on very fast dynamical time scale and the system does not 
have time to explore all of the phase space to find the most 
probable configuration JArad fc Lv ndcn-Bcll 20051 ). 

It was recently found by iLevin et al.1 ( 20081 ) that the 
Lynden-Bell approach, considering the system confined in 
a finite box, is able to quantitatively predict the one particle 
phase space distribution when the out of equilibrium initial 
state is close to the virial requirement, i.e. —1.2 < bo < —0.8, 
where 



60 = 



2K 

Wo 



(1) 



is the initial (i.e., at time t — 0) virial ratio, while Ko and Wo 
are respectively the initial kinetic and potential energy. The 
Lynden-Bell prediction in a confining box is named "cut-off 
Lynden-Bell" and the cut-off is physically justified by the 
realization that the relaxation must be restri cted to a finite 
region of space (jChavanis fc Sommeriall2008l ). Outside this 
range of 60 values the cut-off Lynden-Bell distribution 
is not able t o describe the sta tistical properties of the re- 
sulting QSS (|Levin et al.ll2008l ). When the cut-off is taken 
to infinity the Lynden-Bell distribution is made of a fully 
degenerate Fermi core and particles at infinity, without the 
halo. 

The cut-off Lynden-Bell distribution was found 
to be successful to explain properties of QSS formed 
in one-dimensional gravitating system s, for initia l con- 



ditions near the virial equil i brium (Yamaguchi 20081; 
Joyce fc Worrakitpoonponll2012l;IW orrakitp oonpon fc Joyce! 



20121) . Recently iTeles. Levin fc Pakterl l|2012l ) introduced a 



novel statistical mechanical approach that can avoid some 
of the fundamental assumptions of the Lynden-Bell theory, 
namely ergodicity and phase-space mixing which are gener- 
ally not satisfied for systems with long range forces. 

An interesting attempt to construct a statistical 
mechanics modeling of the vi olent collapse wa s deve l- 
oped in series of pape rs b y IStiavelli fc Bertinl II 19871); 
iBertin fc Trentil (|2003l ); iTrenti. Bertin fc van Albadal 
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IJ2005F ); iTrenti fc Bertinl (|2005l . 120061 ). This provides phys- 
ically motivated distribution functions derived from the 
Boltzmann entropy conserving mass, energy, plus a third 
quantity Q. The problem is, in general, to determine to 
what extent the three quantities are indeed conserved 
during the collapse, i.e. whether the virialized structure 
formed after the collapse have the same number of particles, 
energy and Q of the initial mass distribution. 

Given that the theoretical problem is very difficult, 
one needs to use gravitational N-body simulations as a 
means to perform simple and controlled numerical exper- 
iments. Previous studies of the relaxation of an isolated 
system sta rting with c o ld enough initial c onditions, i.e 
foo 



0. dHenonl Il964l: Ivan Albadal Il982l: lAarseth et all 



19881: iBoilv fc A thanassoula 20061: JPayid fcTheunsl 1 19891 : 



Therms fc Davidlil990l ; I Joyce et aLll2009l) have shown that 
the system undergoes to a large contraction, reaching a min- 
imal size, approximately at td, that scales as r m i„ ~ N^ 1 ' 3 
with the number of points N (at fixed volume V and mass 
density pa — mN/V). This behavior can be explained by 
considering the growth of density perturbations in th e col- 
lapsing phase (|Aarseth et alj 11988 : I Joyce et all 120091 ). By 
neglecting boundary effects, one may treat the problem by 
using the linear approximation of the self-gravitational fluid 
equations in a contracting universe. The minimal radius re- 
sults to be of the order of the unique length scale character- 
izing the system, i.e., the initial average distance between 
nearest neighbors £ « r m i„ oc N -1 ^ 3 . 

It was then shown (jjovce et al.ll2009l ) that a fraction of 
the particles are ejected from the system because during the 
collapse phase they gain enough kinetic energy. The energy 
ejected grows approximately as N 1 ' 3 while the fraction of 
the mass ejected slowly changes with N. The mechanism of 
ejection rises from the interplay of the growth of perturba- 
tions with the finite size of the system. In particular, parti- 
cles lying initially in the outer shells of the system develop a 
net lag of their trajectories compared with their uniform col- 
lapse ones. This lag propagates into the volume during the 
collapse phase and particles in the outer shells gain positive 
energy by scattering through a time dependent potential of 
an already re-expanding central core. The resulting density 
profile of the virialized state is characterized by a power-law 
profile of the type n(r) ~ r~ 4 for r > r c . Interestingly, this 
same profile was found conside ring several different systems 
IStiavelli fc Bertinl J 1984 ll987T l. Note that ejection of mass 
and energy implies that the mass and energy of the virial- 
ized structure are smaller than the total ones, i.e. there is 
no mass and energy conservation in the collapse. 

In this paper we aim of understanding the origin of the 
n(r) ~ r~ 4 density profile, investigating the properties of 
the initial conditions necessary to obtain such a behavior. In 
Sect[2]we briefly review recent studies of the warm and cold 
collapse. The first is defined for the case in which the initial 
virial ratio is close to foo ~ — 1 while for the second close to 
foo w 0. We motivate the physical reasons for such a distinc- 
tion and we present in Sect[3] the results of some N-body 
simulations where we used the same number of particles but 
we have varied bo in the range [—1,0], with uniform space 
and velocity distributions (i.e., water-bag initial conditions). 
We show that there is a clear differences between the struc- 
tures formed when &o < fo§ ft* — 1/2 and foo > &§• We refer 
to these two relaxation processes, respectively, as mild and 



violent: in the latter case the system contracts by a large 
factor, while in the former it approximately maintains its 
original size. In Sect |4] we discuss in detail the case of mild 
relaxation showing that the predictions of the Lynden-Bell 
theory with a cut-off agree well with simulations. Then in 
Sect[5] we show that the main feature of the bo > fo§ case is 
the n(r) ~ r~ 4 density profile, i.e. the formation of a dense 
core and a dilute halo described by such a power-law pro- 
file. In order to explain the origin of this profile we introduce 
a simple and well-motivated physical model in Sect (6] Then 
we discuss (SectO the origin of the critical value fog ~ —1/2. 
Finally we draw our main conclusions in Sect [8] briefly dis- 
cussing the relation with the halo structures observed in 
cosmological N-body simulations. 



2 VIOLENT AND MILD RELAXATION 

As already mentioned, the properties of the QSS resulting 
from the collapse of an isolated self-gravitating, spherical, 
uniform cloud of particles depend on the initial conditions. 
In the literature there have been mostly studied two different 
cases, i.e. with initial virial ratio foo ~ — 1 and foo = 0, that 
we are now going to review in this section. 



2.1 Lynden-Bell theory in a confining box 

Gravitational systems do not reach a time independent 
equilibrium in the thermodynamics sense. Thus the fine- 
grained distribution function of positions f a nd velocities 
v, f(t ,v,r), never reaches a stationary state. iLvnden-Belll 
(|1967l ) developed an approach based on the idea that a 
coarse-grained distribution function f(t,v,r), averaged on 
microscopic length scales, relaxes to a meta-equilibrium dis- 
tribution f(v,r). The statistical properties of such a state, 
differently from the ordinary equilibrium state characterized 
by a Maxwell-Boltzmann distribution, explicitly depend on 
the initial distribution fo(v,f) = f(t — 0, v, f). Lynden-Bell 
argued that the collisionless relaxation should lead to the 
density distribution of levels which is most likely, i.e. the one 
that maximizes the coarse-grained entropy, consistent with 
the conservation of energy, momentum and angular momen- 
tum. 

If the initial distribution is a water-bag, i.e. positions 
are constrained in r £ [0, Rq] and velocities in v G [0, Vb], 



fo(v,f) = rnQ(Ro - r)Q{V - v) 



(2) 



where Q(x) is the Heaviside step function and r/i — 
rji(Ro,Vo) is a constant, th e maximization p rocedure gives 
a Fermi-Dirac distribution (|Levin et al.ll2008l ) 



f(t,v,f) =i] 1 p(v,r) 



'/i 



exp \p(e(v, r) -»)] + ! 



(3) 



where e(v, r) is the mean energy of particles, /3 and /x are 
two Lagrange multipliers required by the conservations of 
energy and the number of particles, 



/ 



d 3 rd 3 f(t,v,r)e(v,r) = eo 



d 3 rd 3 f(t,v,r) = 1 



(4) 
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where eo is the energy per particle of the initial distribution. 
In this context, the incompressibility of the Vlasov dynamics 
plays the same role of the Pau li exclusion principle (see e.g., 
IChavanis fc Sommerial (|2008|)). Then, the density profile is 
simply 



I 



n(r) =N / f(v,r)d 3 v. 



(5) 



In practice, however, what is found is that self- 
gravitating systems usually relax to structures character- 
ized by dense cores surrounded by dilute halos, the distri- 
bution functions of which are quite different from Lynden- 
Bell f(v,f). The failure of the theory was attributed to 
the fact that the violent relaxation occurs on very fast 
dynamical time scale and the system does not have time 
to explore all of the phase space to find the m ost prob- 
able configurations (|Arad fc Lvnden-Belll 120051 ). Numeri- 
cal simulations, starting from out of equilibrium configu- 
ration characterized by an initial virial ratio of bo « —1/2 
also showed that the Lynden-Bell theory, as well as other 
theoretical attempts, are ad odds with numerical results 
JArad fc Johansson! [20051). 

However recently iLevin et al.l (|2008T ) showed that when 
the initial distribution satisfies the virial condition &o ~ — 1 
the system quickly relaxes to a QSS described quantitatively 
by the Lynden-Bell distribution with a cut-off. The cut-off 
originates from the requirement that particles must be con- 
fined in a finite volume of space. The reason for this comes 
from the fact that the possible configurations include those 
in which the mass is distributed throughout space and such a 
configuration dominates the entropy. The Lynden-Bell pre- 
diction in a confining box is referred as "cut-off Lynden- 
Bell". It was then shown that for short enough time scales 
the p recise value of the cut-off is unimportant (| Levin et al.l 
120081 ). The metastable Lynden-Bell distribution persists un- 
til a fraction of the particles evaporates because of two-body 
collisions. 

A similar agreement between the cut-off Lynden-Bell 
distribution and numerical simulations was also found for 
initial conditions close enough to the virial condition, i.e. 
— 1.2 ^ bo ^ —0.8 while outside this range the situation 
drastically changes and the Lynden-Bell distribution is not 
able to describe the statistical properties of the resulting 
QSS. Particularly, this occurs when a fraction of the particles 
can gain enough kinetic energy to be ejected from the system 
in a short time scale, while another part, which remains 
bound, form a dense central core and a dilute halo. This 
latter problem is addressed in the following section. 

It is interesting to note that, when the cut-off of the 
truncated Lynden-Bell distribution is extended to infinity, 
then the distribution function splits into two domains, a 
compact core with zero temperature plus an evaporated frac- 
tion of zero energy particles at infinity. The distribution 
function o f the core is given by that of a fully degenerate 
Fermi gas (|Levin et al.ll2008l ). A detailed comparison of the 
Lynden-Bell theory, including density profiles, velocity and 
energy distributions, with numerical sim ulations in one and 
three spatial dimensions is presented in IWorrakitpoonponl 
J2011I ). 



2.2 Mass and energy ejection 

In this section we briefly summarize the main find- 



ings bvlAarseth et all (ll988l);lBoilv fc Athanassoulal (|2006t ): 



iBoilv et al.l ( 20021 ); Ijovce et al.l ( 20091 ) concerning the col- 
lapse of a cold, uniform and spherical cloud of self- 
gravitating particles. In the idealized limit of an exactly 
uniform spherical distribution different shells do not over- 
lap during the collapse. The radial position r(t) of a test 
particle initially at ro is simply given by the homologous 
rescaling 



r(t) = R(t)r 



(6) 



where the scale factor R(t) may be written in the standard 
parametric form 



*(0 = — (£ + sin(£)) 

7T 



(7) 



and 



TD 



-3-7T 



32Gp ' (8) 

Eqs[5][H] describe the unperturbed spherical collapse model 
(SCM) trajectories. At the time td the system collapses 
into a singularity. In a physical situation the collapse is reg- 
ularized by perturbations which are present in the initial 
conditions at any finite N. At first approximation, one may 
neglect the effect of the boundaries on the evolution of the 
density perturbations, i.e. one can consid er the limit of an in - 
finite (i.e., _Ro — > 00) contracting system (| Joyce et alj|2009h . 
One can then consider the fluid limit and solve the appro- 
priate equations perturbatively as it is usually done in cos- 
mology for an expanding (ra ther than contracting as in this 
case) univers e (IPeeblesll 19801). A m ore detailed approach was 
developed bv lAarseth et al.l (|1988l ) taking explicitly into ac- 
count the system finite size. 

When particles are initially randomly distributed (i.e., 
with Poisson fluctuations) one finds that during the 
collapse t he structure reaches a minimal radius which 
scales as (lAarseth et al. 1988 ; Boilv fc Athanassoulal 120061 ; 



IBoilv et al.ll2002l ; Ijovce et alfcosl ) 



,„ oc Ar~ 1/3 



(9) 



This scaling with N is obtained by simply taking the cri- 
terion that the SCM breaks down when fluctuations at 
a scale of order of the size of the system go non-linear. 
EqO has a very simple interpretation. Neglecting the finite 
size of the system, and given that gravity has no intrinsic 
length scale, on purely dimensional grounds we have that 
Train should be proportional to the only length scale in the 
problem, the mean inter-particle distance £ o c A r ~ 1 ' 3 . Eq[9| 



has b ee n observed in N-body si mulat io ns bvlAar seth et al 
(1 19881); IBoilv fc Athanassoulal (|2006l ); IBoilv et al.l (|2002l ) 
iJovceet al.l(|2009l ) 



It was then noticed by I Joyce et al.l (|2009l ) that, while 
all particles start with a negative energy, after the collapse 
a finite fraction ends up with positive energy which may 
escape from the system. This transfer of energy occurs in 
a very short time around td and depends on N; scaling 
behaviors with the number of particles are manifested by 
the amount of ejected energy and particles. Eq[9l together 
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with some simple approximations which have been tested to 
be valid in the simulations, is the key element to understand 
the observed scaling behaviors. 

The radial density profile of the virialized structure 
formed by bound particle s after the collap se was found to 
have the functional form (|Jovce et al.ll2009f) 



n(r) 



1 + 



(*) 



(10) 



where r c and n c are param eters d e pendin g on TV and 
C = 4 in agreement w i th Irlenonl (I1964T1; Ivan Albadal 
<l 19821); IStiayelli fc Berth! (| 19871 ); iBertin fc Trentil (|2003h : 
IRov fc Perez! 1 20041 ). Simple scaling arguments show that 
r c oc TV -1 ' 3 and n c oc TV 2 . In addition it was also noticed 

max Tq ^^ i min- 

Concerning the mechanism of mass ejection it was found 
that there is a very clear systematic correlation between 
particles initial radial position and ejection, a fact that has 
lead to understand that the physical mechanism of ejection 
indeed arises from the coupling between the evolution of 
pertu rbations and the finite size of the system (jJovce et al.l 
120091 ). Given the importance of such a mechanism for the 
rest of the paper, let us describe it in some details. 

The key to understand the ejection mechanism is to re- 
alize that particles initially lying in the outer boundary lag 
behind the others during the collapse. This lag can be un- 
derstood as follows. Local density fluctuations modify the 
SCM trajectories (i.e., Eqs l6l8|) so that the contraction is 
no more perfectly homologous. In this situation there is an 
asymmetry between the shell at the outer boundary com- 
pared to the ones in the bulk: as particles move around 
there is no compensating inward flux at the boundary for 
the mass which moves out under the effect of perturbations. 
For this reason the density of the outer shell decreases, and 
also the average density in the sphere at the corresponding 
radius, slowing its fall towards the origin. As time goes on 
this asymmetry propagates into the volume and for this rea- 
son particles in the outer shell particles arrive at the center 
of mass on average much later than those in the bulk. 

The mechanism of the gain of energy leading to ejection 
is simply that the outer particles, arriving later on average, 
move through the time dependent decreasing mean field po- 
tential produced by the re-expanding inner mass. It is pos- 
sible to work out a simple estimate for the eje cted energy 
that agrees quite well with the observed scaling l|Jovce et al.l 
120091 ). 

With resp ect to the prediction s of the theoretical model 
introduced by iLvnden- Belli (|1967l ) . it is interesting to note 
that, because of ejection, energy and mass are not conserved 
during the collapse. As discussed in Sect l2.ll this situation 
violates the energy/mass constraints on the final state that 
is assumed in the Lynden-Bell treatment. For this reason, 
it is not surprising that this approach cannot successfully 
explain the statistical properties of the resulting virialized 
structure. 



3 N-BODY SIMULATIONS 

3.1 Initial conditions 

The initial conditions of the simulations are generated as 
follows. We randomly distribute TV particles, of mass m, in a 
sphere of radius Ro with mass density po — 37V/(47ri?o) • rr\j. 
The gravitational potential energy at time t is 



W(t) 



N JV 



Gmiirij 



(11) 



i=i j=i 



where Tij is the distance of the i from the j particle. The 
total kinetic energy is simply 



1 N 



(12) 



i=l 

■th 



where Vi{t) is the velocity of the i particle. The virial ratio 

is 

2K{t) 



b(t) 



W{t) 



(13) 



We generate a series of spherical clouds of particles, 
with TV = 10 4 and with different initial virial ratio bo — 
b(t — 0). We take the velocity components to have a uniform 
probability density function (PDF) in the range [— Vb,Vb], 
and the modulus of the velocity is constrained to be in a 
sphere of radius Vo- The velocity PDF is thus 



g{v) = — 3« 2 for v < Vo 
v o 

and zero otherwise. Such a PDF clearly satisfies 

/■Vo 

g(v)dv = / g(v)dv = 1 . 
Jo 

The initial velocity dispersion is 

rV 



2 (2 3 2 

(v ) = / v g(v)dv = -V 



where we defined 



V,t 



bpGNm 
Ro 



(14) 



(15) 



(16) 



(17) 



To obtain Eq |17l we used that the gravitational po- 
tential energy of a unifor m spherical mass distribution is 
|Binnev fc Tremaind[l994T ) 



Wo 



3 G{mN) 2 
"5 Ro 



(18) 



The initial conditions are thus constrained in a water-bag 
distribution. 



3.2 Code and numerical parameters 

To run N-body simulations we have used the par allel ver- 
sion of the publicly availa ble tree-code GADGET (|Springell 
120051 ; ISpringel et al.ll200ll ). There are various parameters of 
the code that must be tuned in order to have a good accu- 
racy in the time integration: as a control we have used 
both energy and angular momentum conservation, which 



1 Our units are such that po = 1 gr/cm 3 so that tjj 
seconds 



2100 
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are a sensitive mon itoring of the accuracy of the simula- 
tion (|Aarsethl 120031 ). We used a force softening such that 
e/£ = 0.007 where £ w 0.55(47r^/3iV) 1/3 is the initial inter- 
particle distance. Note that the minimal radius r m ;„ of the 
structure in the case of maximum contraction, i.e. when 
&o = 0, is found to b e r m i„ sa £ (see Sect l2.2[l . As discussed 
in I Joyce et al.l (|2009r ). where a number of tests with differ- 
ent values of e were performed, the dynamics of the collapse 
phase and the formation of the QSS remains unchanged as 
long as e < £,r min . 

In addition to the softening length, the accuracy of a 
GADGET simulation is determined by the internal time- 
step accuracy and by the cell-opening accuracy parameter 
of the force calculation We chose the time-step criterion 
of GADGET with r\ = 0.01. In the force calculation we em- 
ployed the new GADGET cell opening criter i on with a high 
force accuracy of qj? = 0.001 (|SpringeJl2005l ; ISpringel et al.l 
l200lf) . 

The behavior of the energy conservation is shown in 
Fig[Tl we have that AE(t)/E «5x 10~ 3 (where E a = 
Wo + Ko is the initial total energjQ) when bo = 0, in the 
range of time we have considered ^ t ^ 4td ; in the other 
cases energy conservation is about ~ 10 -3 . One may note 
that the larger is bo the less accurate is energy conservation 
as the system size gets smaller and particles gain higher ve- 
locities. The latter is the reason for the largest deviation 
in the energy conservation seen for bo = at t « 4td. 
Moreover, the behavior as a function of time of one com- 
ponent (for instance along the x-axis) of the total angular 
momentum shows that it is well conserved during the time 
integration (see inset panel of Fig[TJ . 
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Figure 1. Behavior of the total energy normalized to its initial 
value as a function of time for different values of bo . In the inset 
panel it is shown the behavior of one of the components (i.e., L x ) 
of the total angular momentum as a function of time. 
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3.3 Global behaviors 

The virial ratio as a function of time b{t) shows a different 
behavior depending on b(t = 0) = bo (see FigfJ]). For bo < 
— 1/2, b(t) presents a series of damped oscillations around 
the asymptotic value — 1. Instead, for bo = it presents a 
sharp change of behavior at td ■ In addition, one may note 
that, for t > td, the virial ratio of the fraction of particles 
with negative total energy stabilizes, as expected, around 
bneg ~ — 1, while the virial ratio of all the N system particles 
reaches the an asyiuptotic value that is btot < b neg . 

This behavior is easily explained by considering the 
ejection of a fraction of the particles from the system — i.e., 
for 60 > — l/2 a certain fraction of the particles gain positive 
energy during the collapse. Their kinetic energy is the ori- 
gin of the offset between b to t and b neg . Indeed, the potential 
energy of the particles with positive energy becomes negli- 
gible (i.e., IVKposI <S |Wneg|) because their distance from the 
structure rapidly increases, so that at first approximation 
we have 

btot = -— — « b neg + — < b neg . (19) 

Wtot Wne g 

On the other hand, for bo < —1/2 all particles remain 
bounded to the structure and thus b„ eg (t) = btot(t). 
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Figure 2. Behavior of the virial ratio for all system particles 
(black line) and only for particles with negative total energy (red 
line) as a function of time for different values of the initial virial 
ratio. Upper left panel: bo = 0, upper right panel: bo = —0.3, 
bottom left panel: bo = —0.7 and bottom right pane:l bo = —1. 



2 In the computation of the gravitational potential energy we 
have taken into account the shape of the gadget softened potential 



iSprin gcl 2005; Springel et al.ll2001I) . 
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Figure 3. Fraction of the particles with positive energy function 
of time for different values of bo: for bo < — 1/2 there is no ejection 
of particles. 



Figure 4. Gravitational radius of the structure as function of 
time for different initial virial ratio bo. 



This picture is conformed by Fig[3]that shows the frac- 
tion of particles f P (t) with positive energy as a function of 
time: we find f p (t) > for t > to and bo > —1/2. On 
the other hand, for bo < —1/2 there is no ejection and 

/*(*) = Vt. 

As long as the spherical structure has uniform density 
the gravitational radius 



R B (t) 



3GM 2 



5W(t) 



(20) 



coincides with the physical radius. From the analysis of the 
behavior of R g (t) shown in Fig(?]we may conclude that min- 
imal size of the structure also depends on bo- In particular, 
the minimal size r m i n <C Ro is reached when bo — > 0. while 
for bo < — l/2 the size of the structure is almost unchanged. 

In summary we have found that there is a clear differ- 
ence between the behaviors of the relevant physical quan- 
tities for different initial virial ratio, particularly when the 
bo is smaller or larger than b§ ~ —1/2. In what follows we 
will study the statistical properties of the resulting quasi- 
equilibrium structure: we firstly, in Sect(4] discuss the prob- 
lem of "mild relaxation", i.e. bo < &§, to then pass in Sect[5] 
to the problem of "violent relaxation" for bo > b(j. In Sect[7] 
we will consider the problem of understanding the origin 
the (approximate) value of b§. 



mass) do not change qualitatively the results discussed be- 
low. The density profile is shown in Fig[5] (upper left panel) 
together with the cut-off Lynden-Bell distributiorjj (see 
Sect l2.1[) . which nicely fits the measured behavior. A more 
detailed comparison of the results of N-body simulations 
wi th the predictions o f the Lynden-Bell theory can be found 
in IWorrakitpoonponl (|201lh . where it is discussed that also 
the energy and velocity distributions are in good agreement 
with the theoretical behaviors. The density profile can be 
best-fitted by a function of the type 



n(r) = n c exp(— (r/r c ) v ) 



(21) 



where r\ (a 2. In addition we find that the characteristic 
length scale r c is of the same order Ro, implying that the 
system has not gone through a drastic change of shape and 
size. Rather it is only slightly changed so that particles have 
rearranged their positions and velocities to find a quasi- 
cquilibrium configuration. In Fig[5] (upper right panel) it is 
also shown the behavior of the energy e p of the i th particle 
as a function of its distance from the center. We may note 
that e p < Vi, which corresponds to the fact that all par- 
ticles are bound: note no clear correlation between energy 
and spatial position is detected. 

For an isotropic radial density profile, p(r), one may 
solve, analytically or numerically, the Jeans e quation to 
get the corresponding velo city dispersion, <r 2 (r) (|Hernquistl 
ll99(J : iTremaine et al.lll994T h The Jeans equation is 



4 MILD RELAXATION AND THE 
LYNDEN-BELL PREDICTIONS 

Let us firstly discuss the case bo — —1. Hereafter, we identify 
the center of the structure as the point in which the poten- 
tial is minimum: alternative definitions (i.e. the center of 



1 d(v r (r) p(r)) sVrir) 
■ + a(r) 



p(r) 



dr 



dr 



-n ■ (22) 



3 I thank Yan Levin and Renato Pakter for their data on the 
cut-off Lynden-Bell distribution. 
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Figure 5. Behavior of some statistical quantities for the case 
bo = — 1- Upper left panel: density profile together with the pre- 
diction of the Lyndcn-Bcll theory (LB) and the best fit with 
Eu l21l Upper right panel: energy per particle e p as a function 
of the distance from the center at t = Att). Bottom left panel: 
mean square value of the radial velocity together with the pre- 
diction of the Jeans equation (Eq[24j. Bottom right panel: phase 
space density p/o- 3 (r). 



In the previous equation a(r) 
persion in the radial direction, 

a(r) = 2 



v r (r) is the velocity dis- 



vt(r) 
v r (r) 



(23) 



is the the anisotropy parameter and Vt(r) is the velocity in 

2 2 

the transversal direction. When Vt(r) = v r (r) the velocity 
anisotropy terms are zero and Eq |22l can be rewritten as|j 



Vr{r) = — r^ 

p(r) 
with the boundary condition 



p(y)GM(y) 



dy 



lim v, 



2 

(r) p(r) 







(24) 



(25) 



It is interesting to note that the Jeans equation (Eq |24[l 
is reasonably well satisfied in the time range we consider 
(Fig[5] — bottom left panel): this implies that the stationary 
state is well described by a stationary solution of the Vlasov 
equation, i.e. it is a collisionless stationary state. It should 
be noticed that although the velocity anisotropy (Eq |23l) is 
different from zero (see below) , the perturbation to the Jeans 
equation due to such a term does not sensibly affect the 

2 

agreement between the measured v r (r) and Eq |24l (We will 



4 A more detailed study of the stationary solutions of the Vlasov 
equation should consider the solution of Eq|22| with a non-zero 
anisotropy term (see e.g., lTrenti fc BertirJ 1120061) ). 
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Figure 6. Other statistical quantities for the case 60 = - 1. Upper 
panel: particle energy fluctuations (sec Eq |26t . Bottom left panel: 
energy distribution of all particles at different times. Bottom right 
panel: radial velocity distribution of particles at different times 
together with a the best fir with a Gaussian function. 



come back on this point in Sect[5]) Finally we note that the 
phase-space density p/a 3 (r), where a 2 = (v 2 (r)} is about 
flat, with a sharp decay for r — > Ro- 

A statistical measure of the amount of energy that all 
particles have exchanged can be defined as follows 



<A 2 (i)> 






1 i^j 



(ej(t) - ei(t)) 2 



N(N - 1) 



Z (t)y 



(26) 



where e\,(t) the average energy per particle is defined as 



(e(*)> 



TZxMl 



N 



(27) 



One may see from Fig[6]that (A 2 (i)) oscillates in phase with 
the virial ratio (see FigO and that the amount of energy 
exchanged by all particles is smaller than 10% during the 
whole time range considered. 

The case 6q = —0.7 does not show substantial differ- 
ences with respect to the b — — 1 case (see Figs l7l8|l . The 
prediction of the Jeans equation for the velocity disper- 
sion shows again that the system is well described by the 
collision-less limit (neglecting the term a(r) in Eq |22[) . The 
density profile is still characterized by a constant behavior at 
small scales followed by a sharp decay of the type described 
by Eq[2l] although r c is smaller than for the 60 = — 1 case, 
implying a larger contraction during the collapse phase. Cor- 
respondingly, particle energies, for t > 2td, are larger than 
for the 60 = —1 case, but still e z p < Vi. The exchange 
of energy among particles (Eo 126)1 was more efficient during 
the first oscillation of the system, i.e. for < t < 1.5td, and 
it is then reduced at later times, in agreement with the fact 
that the system is relaxed into a QSS: each particle move in 
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Figure 7. As Fig[5]but for the case 6q = —0.7. 
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Figure 8. As Fig[6]but for the case 6q 



-0.7. 



a time independent potential and the energy of each particle 
is conserved modulo two-body collisions. 



5 VIOLENT RELAXATION AND THE 

FORMATION OF THE POWER-LAW TAIL 
OF THE DENSITY PROFILE 

We now present the main results of N-body simulations for 
the case in which the initial virial ratio is — 1/2 < bo ^ 0. In 
this case during the collapse the size of the system undergoes 
to a large compression and a fraction of the particles gain 
a certain amount of kinetic energy so that they will have 
velocities larger than the escape one. 

In Fig(9] (upper left panel) it is shown the density profile 
at t > td'- one may note that an almost asymptotic behavior 
is reached already at t > Trj. However, at later times the pro- 
file is almost identical but for the fact that the tail extends to 
larger scales. We find that the density profile is well approxi- 
mated by Eq[TD]where r c S3 0.03i? o and £ = 4. Note that the 
density profile has two different regimes: at small scales, i.e. 
r < t c , the structure corresponds to an homogeneous sphere 
with constant density n(r) w n c while at large scales, i.e. 
r > r c it shows a r~ 4 decay. As already mentioned we find 
that r c ~ I — 0.557?o(47r/7V) 1/3 . The minimal radius r min 
reached by the structure during the collapse can be defined 
as the radius, measured from the center of mass, enclosing 
the 90% of the mass. It is found that r m in ~ r c ~ I- 

In the upper right panel of FigO it is shown the dia- 
gram radial distance-total energy only for the bound par- 
ticles. One may note that for r > r c the points follows an 
approximate e p ~ r _1 behavior. This can be easily explained 
by considering that particles at distances r > r c move in a 
constant gravitational potential generated by particles with 
r < r c . In this situation particle velocities should display a 
Keplerian behavior v r ~ r -1 ' 2 so that that e p ~ i; 2 . ~ r _1 . 
This behavior is confirmed by considering (bottom right 
panel of Fig[5} the behavior of the average radial compo- 
nent of the velocity as a function of the radial distance. 
(The average has been performed in radial shells). Indeed, 
we find that 



<«r(r)> 



a 2 (r) 



1 + 



i) 



(28) 



where <r 2 is a constant and r c has been determined from the 
density fit (Eq fT0)l . The behavior of Eq[28]is similar to the 
one of the density in Eq llOl it is constant at small scales and 
it decays for r > r c . As time passes, a few particles reach a 
larger and larger distances from the center of the structure, 
leaving however unchanged the functional behavior of Eq |28l 
Given the behaviors of Ea llOl and Eg 1281 we may fit the 
phase-space density with (see the bottom left panel of FigO 



P( r ) 
a 3 (r) 



1 + 



\ 



1 + 



(29) 



5/2 



Thus we find that the phase-space density is pfo oc r 
for r > r c while it is almost flat at smaller scales. 

It is interesti ng to note, as firs tly shown in the pio 



neering paper by van Albadal ( 1982T) and t hen studied in 



detail by I Trent i. Bert in fc van Albadal (|2005l ), that the QSS 
formed after the collapse is dominated, in the outer regions 
where the density scales as n(r) ~ r~ , by radial orbits. This 
is shown by the behavior of a(r) (Eq[23} as a function of the 
radial distance (see Fig |10[) . On the other hand for the QSS 
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Figure 9. As Fig[5] but for bo = 0. The best fit of the density 
profile with Eq llOl is also shown, together with the fit given by 
Eq |28l of mean square value of the radial velocity and the fit given 
by Eg 1291 for the phase space density. 
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Figure 10. Velocity anisotropy for the case &o = and bo = —1. 
When a(r) > the velocity dispersion for bo = is by dominated 
by its radial component while for bo = —1 hy its transversal one. 



obtained starting from 6o = — 1 the velocity dispersion is by 
dominated by its transversal component. 

The behavior of A 2 (£) (Eq |26[) (see the upper panel of 
Fig lll[) shows that in this case particles exchange a substan- 
tial amount of energy in the time range 0.7i~d < t < 1.2td, 
while before and after the central collapse phase the en- 
ergy per particle is very well conserved. Differently from the 
bo < —1/2 case, during the collapse a fraction of the par- 
ticles change their total energy by a relevant factor so that 
some of the particles may gain enough kinetic energy to es- 
cape from the system. The bo > —1/2 case corresponds to 
an almost instantaneous collapse followed by a rapid relax- 
ation toward a QSS. The particle energy distribution (Fig llll 
bottom left panel) shows that some of the particles have in- 
deed positive energy. Finally the radial velocity distribution 
(Fig llll bottom right panel) is reasonably well fitted by a 
Gaussian function. 

Behaviors similar to the ones shown in FigslMTTI are 
found for the case in which the initial virial ratio is bo = —0.3 
(see Figs ll2|j"3|) . However, due a the non zero initial velocity 
dispersion, the collapse is less peaked in time. The density 
profile is again well approximated by Eq llOl but in this case 
r c — 0.2Ro i.e. it is about ten times larger than for the bo = 
case. Correspondingly we find n c (b a = —0.3) < n c (bo — 0), 
i.e. the structure of the QSS is much less compact. Also the 
behaviors of (v% {r)} and of p(r) / a 3 (r) are well described by 
Eas r28ll29l although with different parameters. Finally A 2 (t) 
shows that there is a smaller exchange of energy during the 
collapse phase than in the bo = case, but still much larger 
than for bo < —1/2. In brief, in this case the collapse is less 
violent and the fraction of particles with positive energy for 




Figure 11. As Fig|6]but for b = 0. 



t > td is greatly reduced with respect to the bo = case 
(see Fig[4j. 
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Figure 12. As Fig|12lbut for the case bg 
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Figure 13. As Fig[TT]but for the case bg = -0.3. 



6 THE ORIGIN OF THE POWER LAW TAIL 
IN THE VIOLENT RELAXATION CASE 

We now introduce a simple physical model with the aim of 
describing the dynamics of bound particles with r > r c when 
bo > bg ~ —1/2. We then show that this model allows us 
to capture the essential ingredients that originate the power 
law r~ 4 tail in the density profile (Eq |10|) . 



6.1 A simple physical model 

We suppose that the bound particles with r > r c for t > To, 
are only subjected to the gravitational field of the core with 
mass 

Mc=AlT L i+{™ry drKi Y ncrlm > (30) 

so that the equation of motion for one of these particles is 
simply 



GM C 



d 2 r _ 
H? ~ r 
We can integrate Eq |31l to get 



1 (dr 

2 \~dl 



GM C 



eo 



where we defined 



eo = 



GM C 
ro 



1 2 

2 Vo 



(31) 



(32) 



(33) 



and ro,vo are respectively the initial position and velocity 
at the initial time tg. 

The initial conditions at io are specified as follows. We 
take the origin of the time at to — to , i.e. the time of max- 
imum collapse of the system. In this situation particles are 
confined in a spherical volume of radius r m i n w r c , the min- 
imal radius of the system during the collapse. As particles 
forming the density power-law tail must be bound, their en- 
ergy is negative, i.e. eo = — e p (ro) > at t — To- We make 
the hypothesis that this energy is conserved at later times. 
This hypothesis is both confirmed by the simulation (see 
below) and compatible with the fact that the system is a 
quasi-stationary equilibrium for t > to- Indeed, in a QSS 
— defined in the mean field limit — each particle moves 
in a time independent potential, and therefore has exactly 
fixed energy. In principle, any change of energy is due to fi- 
nite N effects, which are, however, relevant only on a much 
long time scale than the one considered here. Particles ve- 
locities are assumed to be oriented outwards, an hypothesis 
that agrees with the fact that after the collapse the veloc- 
ity dispersion is by dominated by its radial component (see 

FigED. 

Bound particles may have a maximum velocity such 
i.e. for 



that ep = 



M 
V = 



2GMc 

ro 



(34) 



so that also the velocity is bounded in < vo ^ Vg . By 
defining 



th = 



GMc 
2e 



we can rewrite Eq |32l as 
dr 



-clt 



= V2to— = dr] , 



(35) 



(36) 



v2thT — r" 

where the last equality defines the variable r\. Eq |36l has 
solution in a parametric form 



r(/3) = r H (l + sin(/3)) 



m 



th 

v/2^ 



(/? - cob(/3) - fa + cos(^o)) 



(37) 
(38) 
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Figure 14. Behavior of the distance r(t; vo) reached by a particle 
after a time l/2rrj> from the collapse, as a function of its initial 
velocity v S [0,v^] (from Eqs l37l38t . The behavior of EqflOlis 
also plotted. 



where 



A> 



2ro 



(39) 



Thus for P — /?o we have t = and r(0) = ro. Therefore we 
obtain 



r(0) < r 



2GM C 



W) 2 - M 2 



yp, 



(40) 



where the equality holds for the peaks of the sinusoidal func- 
tion. By inverting Eg 1401 and using Eg 1341 we find 



Vo 



utf 



M\2 



2GM C 



M 
= V 



l-'-l 



(41) 



The behavior of the radial distance rasa function of 
the initial velocity < vo ^ Vq 1 , computed from Eos l37ll38l 
is plotted in FigQ3]for t — 3/2t d . Similarly in Fig[TH]it is 
plotted r(t) for different values of the initial velocity. One 
may note that only the high velocity particles, for which 
Vo ~ Vq, may reach a distance of the order of the initial 
system's radius Ro. As time passes, the particles with the 
highest velocity increases their radial distance to r > Ro. 
At a time of the order of f» 3/2td the structure has already 
reached its (almost) asymptotic shape for r c < r < Ro: while 
at larger scales and at later times, a few particles may arrive 
to larger and larger distances. 

6.2 Shape of the density profile 

Let us now compute, under some simple approximations, 
the density profile resulting from this simple physical model. 
We suppose that all particles have, at the same initial time 
td, the same initial position ro ^ r c . In this approximation, 
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Figure 15. Behavior of the distance r(t) as a function of time 
for different values of its initial velocity vq (from Eqs l37l38t . 



given a certain distribution of initial velocities p(vo), we find 
that the radial density profile, for t > 2td (see Fig |15[) , is 
given by 



^"i^ ^^ 



GM c Np(v ) 



4nr 4 



>/K M ) 2 - 



One may note that from Eo l41l we find that vo 
r > r c , so that in this limit Eg 1421 becomes 



n(r) 



GM C N 

4:Tvr 4 v^ 



P(«rf) , 



(42) 



M r 

V IOr 



(43) 



thus showing the r 4 decay in the best fit of measured n(r) 
(EqflDJ. 

Although in the derivation of Eo 142 1 we made important 
simplifications, we now show that the hypotheses used allow 
us to capture the main elements of the problem. We may 
relax these assumptions by allowing that the initial particle 
positions ro also have a certain PDF /(ro). In this case we 
need to integrate Eos l37l351 numerically as follows: 

• We extract the initial conditions [ro, vo], such that < 
r ^ r c and < do ^ Vq 1 , from the assigned initial position 
and velocity PDFs /(ro) and p(vo)- 

• We fix the time t > td at which we compute the profile. 

• We then find from Eo 1381 the value j3 = /3(t, ro, «o) and 
then from Eo l37l r(t'l = r(£,ro,i>o). 

• We repeat the iteration N times and we can thus con- 
struct numerically, from the resulting distribution of dis- 
tances r(t), the density profile at time t. 

As an example (see Fig |16|) we have assumed a Gaus- 
sian velocity distribution p(vo) with zero mean and variance 
(iiQ / ) 2 /a, where a > 1 is a free parameter, and we have 
considered only particles with i>o > 0. In addition we have 
taken a uniform distance ro distribution such that /(ro) 7^ 
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Figure 16. Behavior of the theoretical radial density n(r) com- 
puted with the Monte-Carlo method described in the text at dif- 
ferent times. The normalization is arbitrary. 
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Figure 17. Behavior of the position at time t (black dots t = 
1/4t£>, red dots t = 3/Atd) as a function of the initial position 
ro, averaged in shells, for different values of the initial virial ratio 
(see caption). 
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only for ro £ [0, r c ]. We have tested that no sensible differ- 
ences are detected as long as a is not large enough to get a 
too small value of p(v ): in particular, when Np(vo ) < 1 a 
particle cannot neither be ejected nor form the r -4 tail. In 
principle, one should also consider the correlations between 
ro and vo and the fact that particles have a distribution 
of initial times: however these complications do not signifi- 
cantly alter the result of Eq |43l 

Finally we note that as time passes, for t > td, the 
power-law tail extends to larger and larger scales (see Fig[§]). 
This is simply explained by the behaviors shown in Figs ll4t - 
1151 the largest distance reached by the particles with highest 
velocity also increases with time, although the precise man- 
ner in which this occurs depends in a very detailed way on 
the value on the properties of p(vo) for vo — > Vq 1 . 



7 THE CRITICAL VALUE OF THE INITIAL 
VIRIAL RATIO 

We now consider the question of what determines the critical 
value of the initial virial ratio for having or not ejection 
(and thus the formation of the n(r) ~ r~ 4 power law tail) 
to be Vq « —1/2. We would like to stress that the precise 
value &o must be a function of N, as collisional and discrete 
effects, although re present perturbati ons, are also present in 
the collapse phase |jovce et al.ll2009u . 

As mentioned in Sect l2.2l the mechanism of energy 
and mass ejection is based on the fact that a fraction of the 
particles, and particularly those that lie at the boundary of 
the system at the initial time, lag behind with respect to 
the others during the collapse, i.e. at £ < td. Particles in 
the bulk collapse approximately satisfying the condition of 



homologous contraction (Eq[6|. The question is whether this 
is also satisfied when bo < 0. 

In Fig |17l we plot, for different values of the initial virial 
ratio bo, particle positions at time t (for t — 1/4td, 3/4td) 
as a function of their positions at time t — 0. We have con- 
sidered an average in shells where these are taken at t — 0; 
we then plot the average value in each shell together with 
the r.m.s. error. One may see that for bo — the two curves 
do not overlap, while this marginally occurs for b — —0.3. 
In this case the homologous contraction is a reasonable ap- 
proximation of the collapse. 

For smaller values of the initial virial ratio, i.e. &o = 
—0.7,-1, there is a substantial overlapping of the curves 
at different times, which means that particles originally be- 
longing to different shells interchange their positions. This 
implies that the collapse cannot be anymore approximated 
as homologous because different shells largely overlap well 
before td- The key mechanism of the growth of the time lag 
is thus eliminated when particles have initially high enough 
velocity dispersion, as different shells cross each other well 
before td- Thus particles from the outer shell arrive at dif- 
ferent time at the center and they do not gain the necessary 
energy to escape from the system. 

Finally it should be noted that the analysis presented 
in this section holds only for water bag initial conditions, 
and that a different initial spatial and velocity distribu- 
tion can lead to very diffe r ent b ehaviors. For instance, 
iTrenti. Bertin fc van Albadal (J2005J) have found that initial 
conditions with equal virial ratio but different spatial distri- 
butions, as for example by generating clumpy distributions, 
lead to si gnificantly less m a ss eje ction. A similar result was 
found by ITrenti fc Bertinl (|2006l ) when particles were ini- 
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tially distributed in a shell rather than in a sphere. Given 
that the precise amount of mass and energy ejection is de- 
termined by the coupling of the growth of perturbation with 
the finite size of the system it is difficult to develop a general 
argument which is valid for different initial conditions and 
a more detailed consideration of each case is needed. 



8 DISCUSSION AND CONCLUSION 

The collapse of an isolated, uniform and spherical cloud of 
massive particles interacting by Newtonian gravity repre- 
sents a paradigmatic example for the formation of quasi 
equilibrium states. It is indeed well known, since the ear- 
liest N-body simulations that, when initial velocities are set 
to zero, this system collapses in a relatively short time scale 
td ~ \/Gpo reaching a configuration which satisfies the 
virial theorem b(t) = 2K/W « — 1. The collective relaxation 
process acting on such a short time scale and the statistical 
properties of the formed quasi equilibrium state have been 
considered in this paper both by performing numerical sim- 
ulations and by an attempt to elaborate a physical model 
able to capture the essential elements of the problem. 

In particular, the initial conditions of simulations are 
generated by randomly placing iV particles with average 
mass density po in a spherical volume and characterized 
by different initial virial ratio ^ bo ^ — 1. Thus, only 
two parameters, i.e., N,bo, define the initial conditions: in 
this paper we hav e varied 6o in the range [—1,0] while in 
Ijovce et al.l (|2009f ) we have considered, for the case bo = 0, 
simulations with different number of points N. 

The system thus collapses under its self-gravity and 
then it forms a virialized structure. This does not repre- 
sent an equilibrium state in the thermodynamics sense. In- 
deed, two body collisions, which have a time scale of about 
i"2 ~ td ln( N)/N will cause a slow e vaporation of parti- 
cles from it (|Binnev fc Tremainelll994h . For this reason the 
state formed after td is called a quasi sta tionary state (QSS) 
(|Dauxois et alJl2003tlCampa et alj|2008h . 

By considering N-body simulations with different initial 
virial ratios, we have identified a critical initial virial ratio 
&o w —1/2 separating the formation of two qualitatively dif- 
ferent kind of QSS. When -1 < &o < -1/2 the collapse 
consists in a series of damped oscillations, the first of which 
one has larger amplitude. The system thus approximately 
maintains its original size. The density profile characteriz- 
ing the virialized QSS is well fitted by the predictions of 
the Lynden-Bell distribution with a cut-off by c onsidering 
the system confined in a box (|Levin et alj|2008r ). This is 
characterized by an abrupt decay of the density at a scale 
r c ~ Rq. The Lynden-Bell predictions strictly holds for a 
confined system: however it was found that the theory does 
not depend sensibly on the cut-off value. This approach is 
thus useful to understand the properties of gentle kind of 
collective relaxation which occurs when the systems is ini- 
tially in a configuration which is close enough to the virial 
equilibrium. 

On the other hand, for —1/2 ^ &o ^ the system size 
undergoes to a large compression and a part of its mass and 
energy is ejected. Finite N fluctuations in the initial spatial 
particle configuration generate density perturbations which 
grow during the collapse. When i>o = such a dynamical 



problem can be treated, when boundaries effects are ne- 
glected (i.e. in the limit Ro — r oo) as the growth of perturba- 
tions in a contracting universe. When fluctuations at a scale 
T m m of order of the size of the system go non-linear the col- 
lapse is stopped I Aarseth et al.l 19881 : Boilv fc Athanassoulal 



l2006t iBoilv et al]|2002l ; Ijovce et alj|2009ft 



During such the collapse some particles gain enough ki- 
netic energy that can be ejected from th e system. The ejec - 
tion mechanism was studied in detail by I Joyce et al.l (| 20090 
where it was shown to be related to a boundary effect. Parti- 
cles initially placed close to the boundaries arrive later than 
the others toward the center, moving, for a short time inter- 
val, in a rapidly varying potential field generated by the par- 
ticles which have already inverted their motion from inwards 
to outwards. In this way they gain some kinetic energy, so 
that some particles have positive energy e p > 0. The density 
profile n(r) of the bound system formed after the collapse is 
characterized by a core, where p(r) ~ const, and by an halo 
in which n(r) ~ r~ 4 . The former behavior can be under- 
stood by considering that the distribution function of the 
core is given by that of a fully degenerate Fermi gas: this 
can be obtained again from the cut-off Lynden-Bell distri- 
bution, by letting the cut-off to extend to infinity so that 
particles can move far away from the center. In this case it 
forms a core-halo s tructures, with a dense core a diluted halo 
(jLevin et alj|2008r ). The Lynden-Bell theory cannot however 
be used to derive the properties of the halo, i.e. that the 
radial density decays as n(r) ~ r~ 4 . 

In order to understand the formation of such a power- 
law tail we have introduced a simple physical model based 
on a few ingredients, namely that that: (i) at the time of 
maximum contraction, i.e. t » td, particles are confined 
in a small phase-space region, (ii) particles energy may be 
close to, or larger than, the escape one and (iii) particles 
forming the power-law tail move in a central and constant 
gravitational potential generated by the mass of the core 
M c at r < r c . With these assumptions a density profile with 
a power-law tail is naturally formed. We conclude that the 
behavior n(r) ~ r~ 4 is the typical density profile that is 
obtained when the initial conditions are cold enough that 
ejection of mass and energy occurs. 

The critical virial ratio 6q separating the two situations 
in which the power-law profile is formed, and mass ejection 
occurs, can be understood by considering that when the ini- 
tial velocity dispersion is large enough the contraction is no 
more homologous. Therefore different shells may overlap be- 
fore the final collapse phase at t m td and the mechanism 
underlying the gain of energy for the outer particles cannot 
be working anymore. 

Finally it is interesting to note that cold dark mat- 

1 



ter halos in cosmological simulations ( Navarro et al 



19971 : lHanser 



1996: 



2004; 



Moore et al.lll988l. l200ll: iNavarro et al. 
Navarro et all 120041 : iMerritt et all l2006f ) display a density 
profile such that n(r) ~ r~ x at small scales and n(r) ~ ?-~ 3 
at large scales: these behaviors are not observed to form 
from the simple initial conditions we have chosen |f| Also 
the phase space density has a different shape, decaying as 
p/a 3 ~ r - 1875 a t ^i scales in cosmological simulations 



5 Although iGraham et al.l ( 120061) found a steeper slope at large 
scales. 
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(jNavarro et alj|2004[ ) while it displays a r~ 5 ' 2 behavior only 
at large enough scales, i.e. r > r c , in the case of struc- 
tures formed from the initial conditions we considered (when 
6§ > —1/2). This difference maybe originated by that the 
fact that cosmological halos are formed from more compli- 
cated initial conditions than the case we considered. How- 
ever, one should also consider that cosmological halos are 
formed in a complex backgrounds so that the hypothesis 
that they are isolated structures maybe not be a valid as- 
sumption. In addition, there is a continuous mass accretion 
so that neither the total mass nor the total energy are con- 
served. A more focused study of these features will be pre- 
sented in a forthcoming work. 

I acknowledge Roberto Capuzzo-Dolcetta, Massimo 
Cencini, Umberto Esposito, Andrea Gabrielli, Michael 
Joyce, Yan Levin and Tirawut Worrakitpoonpon for useful 
discussions and comments. 
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